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Abstract 

A  one-dimensional  mathematic  computational  fluid  dynamics  model  of  a  proton  exchange  membrane  (PEM)  fuel  cell  is  presented  in  this  paper 
to  simulate  the  transient  behavior  of  hydrogen  pressure  in  the  flow  field  during  a  typical  dynamic  process — the  purge  process.  This  model  accounts 
for  the  mechanism  of  pressure  wave  transmission  in  the  channels  by  employing  the  characteristic  line  method.  A  unique  parameter — pressure 
swing,  which  represents  the  top  value  of  pressure  variation  at  certain  point  in  the  channel  during  the  purge  process,  is  brought  up  and  studied  as  well 
as  the  pressure  drop.  The  pressure  distribution  along  the  channel  and  the  pressure  drop  during  the  purge  process  for  different  operating  pressures, 
lengths  of  purge  time,  stoichiometric  ratios  and  current  densities  are  studied.  The  results  indicate  that  the  distributed  pressure,  pressure  drop  and 
pressure  swing  all  increase  with  the  increment  of  operating  pressure.  With  a  high  operating  pressure  a  second-falling  stage  can  be  seen  in  the 
pressure  drop  profile  while  with  a  relatively  low  operating  pressure,  a  homogeneous  distribution  of  pressure  swing  can  be  attained.  A  long  purge 
time  will  provide  enough  time  to  show  the  whole  part  of  the  pressure  drop  curve,  while  only  a  part  of  the  curve  can  be  attained  if  a  short  purge  time 
is  adopted,  but  a  relatively  uniform  distribution  of  pressure  swing  will  show  up  at  the  moment.  Compared  with  the  condition  of  stoichiometric  ratio 
1,  the  pressure  drop  curve  decreases  more  sharply  after  the  top  value  and  the  pressure  swing  displays  a  more  uniform  distribution  when  the  ratio 
is  set  beyond  1.  Different  current  densities  have  no  apparent  influence  on  the  pressure  drop  and  the  pressure  swing  during  this  transient  process. 
All  the  distribution  rules  of  related  parameters  deducted  from  this  study  will  be  helpful  for  optimizing  the  purging  strategies  on  vehicles. 

©  2006  Elsevier  B.V.  All  rights  reserved. 
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1.  Introduction 

PEM  fuel  cells  are  considered  as  one  of  the  most  promising 
alternative  power  devices  for  both  stationary  and  mobile  appli¬ 
cations.  However,  with  regard  to  the  longevity,  the  fuel  cells 
for  vehicle  use  inevitably  cannot  survive  as  a  long  period  as 
for  the  stationary  application.  The  reason  may  lie  in  the  fact 
that  the  fuel  cells  for  vehicle  use  frequently  experience  dynamic 
conditions  which  accelerate  the  aging  process  of  the  fuel  cells 
[1].  Pressure  wave  transmission  generated  in  the  flow  channels 
during  the  dynamic  processes  may  lead  to  mechanical  vibra¬ 
tions  of  the  membrane  electrode  assembly  (MEA)  and  accelerate 
its  damaging.  [2].  Early  dynamic  model  by  Amphlett  et  al. 
[3]  based  on  coupling  the  steady  state  electrochemical  model 
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with  an  unsteady  thermal  model  predicted  the  transient  elec¬ 
trical  response  of  start-up,  load-step  and  slow-down  operating 
conditions.  Most  of  the  following  dynamic  models  [1,4-6]  in 
the  past  were  also  set  up  to  simulate  the  transient  behaviors  of 
current-voltage  relation  polarization  curve,  or  develop  model¬ 
ing  methods  and  experimental  techniques  for  the  study  of  current 
distribution  under  transient  conditions  [7,8],  whereas  few  cal¬ 
culations  have  concentrated  on  the  dynamics  of  pressure  and 
flow  conditions  or  specially  modeling  the  transient  response  of 
pressure  distribution  in  the  flow  channels.  Though  there  had 
been  many  excellent  models  through  CFD  approach  to  simulate 
the  flow  and  the  along-channel  pressure  drop  [9,10],  most  of 
them  were  limited  to  the  pressure  characteristics  under  steady 
state.  Pathapati  et  al.  [11]  proposed  a  transient  mathematical 
model  incorporating  the  effects  of  dynamics  of  flow  and  pres¬ 
sure  in  both  anode  and  cathode  side  and  charge  double  layer  to 
simulate  the  dynamic  response  of  pressure  as  well  as  voltage  out¬ 
put.  However,  the  dynamic  properties  of  pressure  presented  by 
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them  only  characterize  the  average  pressure  in  channels,  while 
the  pressure  wave  at  certain  point  or  non-uniformity  of  pres¬ 
sure  distribution  could  not  be  caught  since  uniform  pressure  in 
the  flow  field  has  been  assumed.  Um  et  al.  [12]  developed  a 
multidimensional  model  through  CFD  approach  to  simulate  the 
transient  details  of  electrochemical  kinetic,  current  distribution, 
hydrodynamics  and  multicomponent  transport  in  PEMFC.  They 
studied  the  non-uniformity  in  the  flow  field  during  transient  pro¬ 
cesses,  but  as  mentioned  above,  this  study  also  focused  on  the 
dynamic  behaviors  of  polarization  relation  and  local  current  dis¬ 
tribution.  Although  the  along-channel  distribution  of  reactant 
and  water  vapor  fraction  had  been  presented  simultaneously, 
the  concentration  of  this  study  was  still  on  the  basis  of  per¬ 
formance  output,  other  than  the  potential  damages  from  the 
pressure  dynamics  in  channels.  Since  the  membrane  electrode 
assembly  (MEA)  of  fuel  cell  suffers  from  frequent  pressure  fluc¬ 
tuations  generating  in  the  channels  during  dynamic  processes, 
a  comprehensive  understanding  about  the  phenomena  and  reg¬ 
ularity  of  pressure  distribution  is  indispensable  to  improve  the 
controlling  strategies  for  various  transient  processes.  Another 
important  point  that  should  be  mentioned  is  that  numerous 
dynamic  simulations  have  been  proposed  while  experimental 
data  about  the  dynamic  processes  are  scarce  in  the  present. 
This  work  is  to  pursuit  the  dynamic  characteristic  of  pressure 
in  the  flow  field  during  a  typical  dynamic  process — the  purge 
process  based  on  one-dimensional  non-steady  fluid  simulation 
and  experiments. 

2.  Model  development 

The  PEM  fuel  cell  model  presented  here  is  a  one-dimensional, 
transient,  computational  model  providing  a  detailed  description 
of  the  following  phenomena  during  the  purge  process  on  differ¬ 
ent  operating  conditions: 

1.  transient  distribution  of  pressure  along  the  channel  in  the 

anode; 

2.  pressure  drop  of  the  gas  flow  in  the  anode. 

The  equations  governing  the  purge  process  include  the  mass 
and  momentum  conversion  equation  governing  fluid  flow,  pro¬ 
cedure  equations  governing  the  assumed  isentropic  process  and 
the  ideal  gas  state  equation. 

The  governing  equations  and  appropriate  boundary  condi¬ 
tions  were  implemented  and  solved  in  MATLAB  environment 
based  on  the  characteristic  line  method  which  has  been  classi¬ 
fied  as  an  excellent  solution  for  the  one-dimensional  unsteady 
flow  condition  solving  [13]. 

2.7.  Assumptions 

A  complete  computational  fluid  dynamics  model  of  the  flow 
field  in  fuel  cell  will  involve  many  complex  factors  and  related 
equations,  which  make  the  model  complicated.  Therefore,  some 
assumptions  should  be  made  to  simplify  the  model. 

From  a  typical  profile  of  pressure  curve  during  the  purge 
process  shown  in  Fig.  1,  the  pressure  variation  can  be  divided 


Fig.  1.  A  typical  pressure  curve  of  the  purge  process. 


into  three  stages:  trailing  edge,  wave  hollow  and  rising  edge. 
According  to  numerous  wave  forms  collected  from  the  tests, 
consistent  with  the  curves  in  Fig.  1,  it  can  be  concluded  that  the 
rising  stage  has  to  take  relatively  long  time  to  recovery  the  initial 
pressure,  compared  with  the  other  two  stages.  Since  the  recovery 
process  proceeds  smoothly  and  no  fluctuations  of  pressure  has 
been  observed,  the  rising  stage  is  ignored  in  this  simulation  to 
save  the  computational  resource  and  improve  the  efficiency  of 
calculation. 

The  modeled  period  including  the  trailing  edge  and  the  wave 
hollow  lasts  a  very  short  time,  usually  less  than  2  s,  so  the  heat 
transfer  between  the  gas  flow  and  outside  during  this  period  can 
be  ignored  and  an  isentropic  flow  condition  is  assumed  here. 
Thus,  the  energy  equation  can  be  neglected  in  the  model.  The 
adiabatic  index  k  is  assumed  to  be  a  constant  and  k  is  only  related 
to  the  diatomic  gas  H2. 

The  effect  of  the  channel  bends  on  the  pressure  drop  is 
ignored  here  and  the  frictional  is  taken  as  the  total  pressure 
drop  in  the  channel.  The  discrete  pressure  drop  due  to  the 
bends  in  the  channel  is  much  less  than  the  frictional  pressure 
drop  due  to  the  viscosity  of  gas  flow,  and  a  calculated  fric¬ 
tional  pressure  drop  can  be  very  near  to  the  experimental  data 
[14]. 

In  the  channels,  water  is  assumed  to  exist  in  the  vapor  phase. 
This  assumption  is  made  to  eliminate  the  two-phase  flow  condi¬ 
tion  in  the  model.  While  the  liquid  water  may  exist  in  practice, 
this  is  beyond  the  scope  of  this  work. 

Other  assumptions  used  in  this  model  are: 

1.  The  hydrogen  gas  is  assumed  to  be  fully  saturated  by 
water  vapor  and  the  mixture  of  hydrogen  gas  and  water 
vapor  behaves  as  compressible  ideal  gas  in  the  anode 
channels. 

2.  The  flow  is  assumed  laminar  due  to  the  low  characteristic 
Reynolds  number.  The  entrance  and  exit  effect  are  neglected 
due  to  the  length  of  channel. 

3.  Similar  flow  conditions  are  assumed  in  all  of  the  channels. 
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Fig.  2.  Geometry  of  serpentine  flow  channel  fuel  cell.  The  red  single  channel 
will  be  the  computational  domain.  (For  interpretation  of  the  references  to  color 
in  this  figure  legend,  the  reader  is  referred  to  the  web  version  of  the  article.) 

2.2.  Computational  domain 

A  computational  model  of  all  the  channels  in  the  flow  field 
requires  excessively  large  computing  resource,  and  since  the 
similar  flow  conditions  among  all  the  channels  are  assumed,  the 
calculation  can  be  based  on  one  serpentine  channel  shown  in 
Fig.  2,  which  is  shown  by  the  red  lines. 

The  computational  mesh  named  position  diagram  [13]  is 
shown  in  Fig.  3.  The  coordinates  in  the  x-t  plane  represent  the 
distance  from  the  inlet  along  the  channel  and  the  time,  respec¬ 
tively.  The  steps  for  calculation  are  determined  differently  for  the 
two  coordinates:  the  time  step  depends  on  the  length  of  calcu¬ 
lated  time,  or  rather  the  purge  time  and  the  accuracy  needed 
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for  the  calculation;  while  the  distance  step  should  be  deter¬ 
mined  by  the  stability  condition  of  the  characteristic  line  method 
[13].  Here,  the  steps  chosen  for  all  the  calculations  are  0.1  s  and 
0.1545  m  and  the  calculations  have  all  been  obtained  on  a  PC 
with  2.4  GHz  Pentium  IV. 


2.3.  Governing  equations 


In  the  fuel  cell  channel,  the  one-dimensional  non- steady  flow 
is  governed  by  the  continuity  equation: 


dp  d(pu)  _  Q 
3 1  dx 

and  momentum  equation: 


3  u  3  u  1  dp 

—  T  u  —  ~b - =  0 

3 1  dx  p  dx 


(1) 

(2) 


where  p,  p  and  u  are  the  density  of  the  gas  mixture,  gas  pressure 
and  flow  rate  in  the  channel,  respectively. 

Since  the  isentropic  flow  is  assumed  in  this  model,  energy 
equation  can  be  neglected,  but  the  whole  process  modeling  here 
will  be  governed  by  the  procedure  equation  due  to  this  assump¬ 
tion: 

P 

—r  =  constant  (3) 

pk 


The  state  of  all  the  gaseous  species  are  governed  by  the  ideal 
gas  state  equation: 


p  =  pRT 


(4) 


where  R  is  the  universal  gas  constant  and  T  is  the  temperature 
of  the  gas  in  the  channel. 

The  equations  above  make  a  set  of  fundamental  equations  of 
the  model.  With  the  sonic  equation:  a  =  ^dp/dp  =  +J k(p/p ) 
[13]  and  Eq.  (3):  p/pk  =  constant,  Eq.  (1)  becomes 


da  da  k  — 

—  T  u  —  T - 

dt  dx  2 


1  3  u 
-a —  =  0 
dx 


(6) 


and  Eq.  (2)  becomes 


du  du  2a  da 

—  T  u  —  T - 

dt  dx  k  —  l  dx 


(7) 


To  invert  Eqs.  (6)  and  (7)  into  ordinary  differential  equations, 
two  combinations  have  been  made:  Eq.  (6)  d=  (k—  l)/2  Eq.  (7), 
then  the  following  equations  are  attained: 


da  da 

— -  +  (u  +  ot)  — 
dt  dx 


+ 


k  -  1 


da 

Yt 


+  (u  —  a) 


da 

dx 


k-  1 
2 


du  du 

— -  +  (u  +  ot)  — 
dt  dx 


du  du 

iF+(“-o)S 


=  0 

=  0 


(8) 

(9) 


Eqs.  (8)  and  (9)  are  obtained  from  the  set  of  fundamental  gov¬ 
erning  equations,  therefore  the  solutions  of  these  two  equations 
equal  to  the  ones  of  the  set  of  fundamental  equations. 

For  Eqs.  (8)  and  (9),  if 

dx 

—  =  u  zb  a 
dt 


Fig.  3.  Rectangular  grid  on  the  x-t  plane  for  calculation. 


(10) 
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then  Eqs.  (8)  and  (9)  become 


da  —  1  du 

dt  2  dt 


(ID 


Eqs.  (10)  and  (11)  constitute  the  ordinary  differential  equations 
of  the  model.  For  the  calculation  with  the  characteristic  line 
method,  non-dimensional  parameters  are  generally  used,  as  the 
following  defines  [13]: 


A 

Z 


a 
ur 
ur  t 
Ir 


u 

U=—, 

UR 


(12) 


where  aR,  pR  and  /r  are  the  reference  sonic  velocity,  the  ref¬ 
erence  pressure  and  the  reference  length  and  A,  U,  X  and  Z 
denote  the  non-dimensional  pressure,  the  non-dimensional  flow 
rate,  the  non-dimensional  length  and  the  non-dimensional  time. 
Therefore,  the  ordinary  differential  equations  of  the  model  with 
non-dimensional  parameters  become 


dX 

—  =  U±A 
dZ 


d  A  k-ldU 

-  zb - =  0 

dZ  2  dZ 


(13) 

(14) 


2.4.  Boundary  condition  and  model  parameters 


but  the  non-dimensional  pressure  and  non-dimensional  flow  rate 
at  the  inlet  will  change  dramatically  with  the  pressure  impulse 
at  the  outlet  during  the  purge  process. 

The  initial  value  of  A(Z=0)  is  determined  by  the 
equations:  As  =  (ps/pRp~mk,  Ain  =  (pm/pR)(k~l)/2k,  Aout  = 
(Pout/PR)(k~mk  and  pout  = pin  —  A p.  Here,  ps  denotes  the  pres¬ 
sure  in  the  gas  store  and  the  initial  value  of  ps  varies  with  different 
cases:  1.50, 1.75  and  2.00  bar  and  A p  is  the  pressure  drop  of  the 
whole  channel  in  steady  state,  which  can  be  calculated  exactly 
based  on  the  equation  provided  in  ref.  [14].  In  this  work,  the  ini¬ 
tial  values  of  pressure  drop  of  the  purge  process  were  obtained 
from  the  corresponding  tests  in  steady  state.  As  assumed  above, 
the  frictional  pressure  drop  is  taken  as  the  total  pressure  drop, 
thus  the  initial  values  of  pressure  at  other  calculated  positions 
along  the  channel  can  be  obtained  with  the  liner  interpolation 
method. 

The  initial  value  of  hydrogen  gas  consumption  is  given  by 

H2  usage  =  — —  =  5.182  x  10"6  x  /  (17) 

2  x  F 


where  F  and  I  denote  the  Faraday  constant  and  the  current  of  the 
fuel  cell. 

Since  the  saturated  hydrogen  gas  is  assumed  in  the  channel, 
the  molar  fraction  of  water  vapor  is  given  by  ref.  [15] 


(IB) 


From  the  items:  A  =  (a/aR)  =  (p/pR)^k~^/2k  and  U  =  u/a r  in 
Eq.  (12),  the  boundary  conditions  of  this  model  depend  on  the 
pressure  and  the  velocity  at  the  inlet  and  outlet  of  the  channel. 

At  the  outlet,  pressure  boundary  condition  is  applied  and  the 
impulse  of  pressure  during  the  purge  process  is  described  as  the 
following: 


A  2k/(k—l) 
^oui  PR 

2k/(k—l) 
Lout 


x2k/(k-\) 


*0 

,  2k/(k—l) 

LR 


Pr  bvt ,  0  <  t  <  tp 


t  >  tr> 


bV  = 


x2k/(k-\) 

20 


—  A 


2k/(k—l) 


(15) 


-PR 


where  Ao  and  Ar  denote  the  initial  and  the  ending  values  of 
the  non-dimensional  pressure  at  the  outlet,  Aout  stands  for  the 
value  of  the  non-dimensional  pressure  at  the  outlet  and  bp  and 
tp  represent  the  slope  coefficient  of  pressure  decreasing  and  the 
time  of  trailing  edge  of  the  purge  process,  respectively. 

At  the  inlet,  a  quasi-steady  flow  condition  is  assumed  here 
because  when  the  hydrogen  gas  is  inducted  from  a  gas  source 
into  the  channel,  the  variations  of  the  parameters  change  more 
dramatically  on  the  v-axis  than  on  the  £-axis.  Therefore,  the  con¬ 
dition  of  parameters  at  the  inlet  in  this  transient  process  can  be 
described  by  the  energy  equation  of  steady  flow  condition 


Al  ~  Ain  + 


-uz 


(16) 


where  As  and  A[n  denote  the  values  of  the  non-dimensional  pres¬ 
sures  in  the  gas  store  and  at  the  inlet,  respectively.  For  every  case 
modeled,  the  non-dimensional  pressure  in  the  store  is  a  constant, 


where  the  subscript  “w”  denotes  the  water  vapor  and  “sat”  means 
saturation  value,  p™1  depends  on  the  inlet  gas  temperature  Tm. 
According  to  the  ideal  gas  state  equation:  V  =  ( nRT)/p ,  the  vol¬ 
ume  flux  of  the  humidified  hydrogen  gas  can  be  given  by 

_  H2  usage  X  fl7jn  _  5.182  x  1( T6//?7jn 
p(  l-xw)  p(  l-xw) 

where  the  subscript  “in”  denotes  the  parameters  at  the  inlet. 

Therefore,  the  initial  value  of  non-dimensional  flow  rate 
UiZ-  0)  can  be  given  by 


u  _  Gin 

aR  dh\/kRTR 


(21) 


where  7r  stands  for  the  reference  temperature  and  d  and  h  rep¬ 
resent  the  geometric  width  and  depth  of  the  channel. 

The  parameters  used  for  the  base  case  simulations  presented 
here  are  shown  in  Table  1 . 


2.5.  Supplemental  conditions  in  algorithm 

The  ordinary  differential  equations  and  appropriate  boundary 
conditions  above  can  be  solved  based  on  the  algorithm  of  the 
characteristic  line  method  to  obtain  the  distribution  of  pressure 
in  the  channel.  However,  in  this  work,  the  compensating  effects 
of  two  regulators — the  relieve  valve  unit  and  the  flux  controller 
on  the  transients  of  pressure  and  velocity  must  be  considered  in 
the  calculation. 

The  relieve  valve  unit  sets  the  value  of  induction  pressure 
and  regulates  the  pressure  when  there  are  some  fluctuations  of 
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Table  1 

Geometrical,  operational  parameters  for  the  cases 


Property 

Value 

Channel  length,  L  (m) 

0.927 

Channel  width,  d  (m) 

0.002 

Channel  height,  h  (m) 

0.0008 

Area  of  the  active  area  of  the  ME  A,  S  (m2) 

154  x 10“4 

The  initial  value  of  temperature,  T  (K) 

323 

Inlet  gas  temperature,  Tm  (K) 

323 

Pressure  in  the  gas  store,  ps  (bar) 

1.50,  1.75  and  2.00 

The  initial  value  of  pressure  drop,  A p  (bar) 

0.004  and  -0.012 

Reference  length,  /r  (m) 

0.927 

Reference  temperature,  7r  (K) 

298 

Reference  pressure,  pr  (bar) 

1.01 

Current,  1  (A) 

92.4,  123.2  and  154 

Hydrogen  stoichiometric  ration,  Za 

1,  1.1,  1.2 

Adiabatic  index,  k 

1.4 

The  universal  gas  constant,  R  (kJ  mol-1  K-1) 

8.314 

Trailing  edge  time  of  the  purge  process,  tp  (s) 

0.07-0.15 

the  pressure  behind  the  valve.  The  flux  meter  works  as  a  flux 
controller  when  the  real  flux  of  gas  flow  exceeds  the  value  set 
previously  and  otherwise  it  is  used  just  as  a  flux  meter.  Therefore, 
when  the  pressure  suddenly  decreases  and  the  flux  exceeds  the 
value  set  beforehand  during  the  purge  process,  the  relieve  valve 
unit  and  the  flux  meter  will  regulate  the  parameters  as  feedback 
controllers.  Since  the  technical  data  with  the  relieve  valve  unit 
and  the  flux  meter  cannot  exactly  describe  their  compensating 
effects  in  this  transient  process,  two  iterative  equations  are  intro¬ 
duced  in  the  calculation  to  simulate  the  contributions  of  these 
two  regulators  in  this  transient  process. 

Px+\  =  p'r+\  +  (Pt  -  P'r+0  X  Ci; 

Mr+ 1  =  u'T+ 1  +  (Mr  -  m'+1)  X  C2.  (22) 

where pr+\  and  ur+\  are  the  pressure  and  the  flow  rate  value  cou¬ 
pling  the  compensating  effects  of  the  regulators,  while  /?'+1  and 
w'+1  are  the  ones  without  coupling  these  compensating  effects, 


the  subscript  “r  +1”  stands  for  the  sequence  of  the  time  step  and 
ci,  C2  are  the  compensating  factors  dependent  on  the  regulators 
themselves.  With  Eq.  (12),  the  two  iterative  equations  become 

/2k/(k—\)  _  A,2k/(k-l)  ( A2k/(k-l)  _  w2*/(*-l)\ 

Ar+1  ~  r+1  r+1  J  X  cl’ 

UT+ 1  =  u'+l  +  (UT  -  U'+1)  X  c2.  (23) 

These  equations  supplement  the  compensating  effects  of  the 
two  regulators  on  the  transient  behaviors  of  pressure  and  flow 
rate  in  the  calculation,  which  improve  the  accuracy  of  the  mod¬ 
eling.  To  simulate  the  characteristics  of  these  two  regulators 
accurately,  numerous  tests  on  different  operating  conditions  are 
necessary  to  decide  the  corresponding  factors  before  the  calcu¬ 
lation. 

3.  Results  and  discussions 

3.1.  Experimental  validation 

To  validate  the  numerical  model  presented  above,  three  com¬ 
parisons  with  different  operating  pressures,  lengths  of  purge  time 
and  current  densities  were  made  between  simulation  results  and 
experimental  data.  The  experiments  were  conducted  on  a  test 
bench  shown  in  Fig.  4. 

As  Fig.  4  shows,  the  pressure  and  flux  of  inlet  hydrogen  gas 
were  controlled  and  stabilized  by  relieve  valve  unit  and  flux 
controller,  respectively.  The  temperature  and  pressure  of  inlet 
and  outlet  gas  were  both  recorded  by  the  collector  with  50  kHz 
frequency.  In  the  terminal  of  the  flow  pipe,  two  outlet  valves 
were  set  and  the  transient  purge  process  can  be  realized  using 
electromagnetic  valve. 

The  pressure  measurements  for  the  purge  process  were  made 
under  the  condition  of  50  °C  cell  temperature.  The  hand  operated 
valve  was  closed  before  the  purge  process,  thus  the  anode  stoi¬ 
chiometric  ratio  was  1 .  This  condition  was  chosen  close  to  the 


„  „  Check  valve 

Pressure  Relieve  valve  Flux 


Hydrogen 
gas  outlet 


Fig.  4.  PEM  fuel  cell  testing  bench. 
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one  under  which  water  accumulations  in  the  channels  easily  took 
place,  followed  by  purge  processes  in  the  real  operating  condi¬ 
tions.  The  reactants  were  fully  humidified  by  setting  a  higher 
dew  point  temperature  than  that  of  the  cells.  It  should  be  noted 
that  significant  water  accumulations  in  the  flow  field  must  be 
avoided  in  the  tests  in  which  two-phase  flow  should  be  consid¬ 
ered.  As  the  assumption  mentions  above,  that  condition  has  been 
out  of  the  scale  of  this  study. 

The  outlet  pressure  profiles  were  attained  through  the  tests 
and  then  adapted  as  the  input  conditions  of  the  simulations, 
accordingly  the  inlet  pressures  can  be  obtained  from  both  sim¬ 
ulations  and  tests.  Figs.  5-7  show  the  comparisons  between  the 
simulation  results  and  the  experimental  data  at  three  different 
operating  pressures,  lengths  of  purge  time  and  current  densities, 
respectively.  Good  agreements  can  be  seen  in  the  comparative 
profiles  of  inlet  pressure.  Inadequate  agreement  for  the  compar¬ 
ative  profiles  of  outlet  pressure  is  attributed  to  the  fact  that  all  the 
outlet  pressure  profiles  were  simplified  into  liner  trailing  edge  in 
simulation,  while  the  real  pressure  at  the  outlet  might  fluctuate 
with  the  pressure  in  pipes  slightly.  But  the  trend  can  be  still  cap¬ 
tured  by  the  present  model.  Pressure  non-uniformity  during  this 
transient  process  can  be  displayed  by  these  simulation  results 
that  are  unattainable  from  experiments. 

3.2.  Pressure  distribution  during  the  purge  process 

A  detailed  distribution  of  local  pressure  in  the  channel  is 
shown  in  Fig.  8.  As  expected,  a  nearly  uniform  decreasing 
trend  of  distributed  pressure  would  result  when  a  low  operat¬ 
ing  pressure  like  1.25  bar  is  adopted,  except  that  some  short 
stable  periods  of  pressure  take  place  at  the  points  (X/L  =  1/2,  2/3 
and  5/6)  in  the  channel  immediately  after  the  trailing  edge.  These 
may  be  caused  by  the  relatively  short  trailing  edge  in  this  experi¬ 
ment,  thus  they  are  the  reflections  of  pressure  wave  superposition 
which  actually  slow  up  the  decreasing  process.  Obviously,  it  can 
also  be  indicated  from  the  figure  that  the  points  in  the  down¬ 
stream  or  close  to  the  exit  of  channel  will  be  influenced  by  the 
trailing  edge  of  pressure  more  noticeably  and  the  stable  pressure 
can  be  attained  more  quickly  than  the  ones  far  from  the  channel 
exit.  Here,  one  unique  parameter — pressure  swing  is  defined 
as  the  highest  amplitude  of  pressure  variation  at  certain  point 
during  the  whole  process  modeled,  and  actually  it  represents 
the  mechanical  attack  of  gas  flow  to  the  membrane  electrode 
assembly  in  this  transient  processes.  The  pressure  swing  versus 
fractional  location  along  the  channel  flow  path  is  given  in  Fig.  9. 
From  Fig.  9,  a  nearly  homogeneous  pressure  swing  distribution 
can  be  seen.  This  is  due  to  the  fact  that  both  the  local  highest 
and  the  lowest  pressure,  actually  the  initial  and  the  stable  value, 
along  the  flow  path  decrease  slightly  during  the  process  when  a 
relatively  low  operating  pressure  is  adopted. 

3.3.  Pressure  difference  effects 

Pressure  distribution  during  the  purge  process  is  influenced 
greatly  by  the  operating  pressure  since  the  ending  pressures  of 
the  processes  for  all  the  experiments  and  simulations  are  set  close 
to  the  standard  atmosphere  pressure.  Fig.  10  shows  the  plots  of 


Fig.  5.  Comparison  of  experimental  and  model  pressure  curves  during  the  purge 
process  for  three  different  operating  pressures:  (a)  1.50  bar,  (b)  1.75  bar  and  (c) 
2.00  bar.  7.a  =  1  at  1.0  A  cm-2;  purge  time,  0.7  s;  operating  pressure,  2.00  bar. 
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Time  (s) 

Fig.  6.  Comparison  of  experimental  and  model  pressure  curves  during  the  purge 
process  for  three  different  lengths  of  purge  time:  1.0,  0.68  and  0.36  s  (from  left 
to  right).  Xa  =  1  at  1.0  A  cm-2;  operating  pressure,  2.00  bar. 

pressure  behaviors  at  different  positions  in  the  flow  field  during 
the  purge  process  for  three  different  operating  pressures.  With 
the  increment  of  operating  pressure  from  1.50  to  2.00  bar,  all  the 
distributed  pressure  curves  are  enhanced  to  new  levels  by  steps 
except  the  outlet  pressure  profile,  which  means  that  the  pressure 
value  at  the  same  instant  for  every  preset  point  will  be  improved 
due  to  the  increased  operating  pressure.  As  can  be  seen  from  the 
figures,  when  the  operating  pressure  exceeds  1 .75  bar,  in  Fig.  10b 
and  c,  a  second-falling  pressure  profile  endings  are  added  to  the 
curves.  The  pressure  variations  at  the  points  near  the  inlet  of 
flow  field  experience  another  decrease  process  after  the  appar¬ 
ently  stable  period  of  pressure  with  a  much  slower  decreasing 
rate,  however  this  kind  of  decrease  in  the  end  of  the  purge  process 
cannot  be  observed  while  the  operating  pressure  is  relatively  low, 
at  1.50  bar.  The  apparently  stable  stage  of  pressure  decreasing 
is  probably  due  to  a  high  pressure  ratio  between  inlet  and  out¬ 


Fig.  7.  Comparison  of  experimental  and  model  pressure  curves  during  the  purge 
process  for  three  different  current  densities:  /=0.6,  0.8  and  1.0  A  cm-2  (from 
left  to  right).  A.a  =  1;  purge  time,  0.64  s;  operating  pressure,  2.00  bar. 


Time  (s) 

Fig.  8.  Pressure  distributed  curves  during  the  purge  process  for  a  low  operating 
pressure  1.25  bar.  Xa  =  1  at  1.0  A  cm-2;  purge  time,  0.7  s. 

let,  which  leads  to  more  pressure  wave  superposition  with  high 
values  and  consequently  slows  the  decreasing  process.  With  the 
weakening  of  the  pressure  wave  amplitude,  the  second-falling 
stage  shows  up.  During  the  terminal  stage  of  the  purge  process, 
or  rather  close  to  the  ending,  the  pressure  drop  between  every 
two  preset  positions  decreases  along  the  channel  for  the  rela¬ 
tively  low  operating  pressure  (1.50  bar),  but  the  trend  changes 
since  the  second-falling  stages  show  up  for  the  relatively  high 
operating  pressures.  In  Fig.  10b,  the  pressure  drops  between  two 
adjacent  points  upstream  in  the  end  of  process  is  mostly  equal 
to  those  downstream.  When  the  operating  pressure  rises  up  to 
2.00  bar  in  Fig.  10c,  the  pressure  drops  between  two  adjacent 
points  upstream  are  much  higher  than  those  downstream.  As  a 
conclusion,  in  the  ending  of  the  purge  process,  the  whole  pres¬ 
sure  drop  will  be  contributed  more  by  the  downstream  field  of 
the  channel  with  the  operating  pressure  increasing,  other  than 
equally  contributed  by  every  area  in  the  flow  field. 


Fig.  9.  Pressure  swing  vs.  fractional  location  along  the  anode  channel  during 
the  purge  process  for  a  low  operating  pressure  1.25  bar.  Xa=  1  at  1.0  A  cm-2; 
purge  time,  0.7  s. 
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Fig.  10.  Pressure  distributed  curves  during  the  purge  process  for  three  differ¬ 
ent  operating  pressures:  (a)  1.50  bar,  (b)  1.75  bar  and  (c)  2.00  bar.  Aa  =  l  at 
1.0  A  cm-2;  purge  time,  0.7  s. 


Fig.  11.  Pressure  drop  curves  during  the  purge  process  for  three  different  oper¬ 
ating  pressures:  1.50,  1.75  and  2.00  bar.  Aa  =  1  at  1.0  A  cm-2;  purge  time,  0.7  s. 


Fig.  11  gives  the  comparison  of  pressure  drop  curves  for 
three  different  operating  pressures.  It  can  be  seen  that  the  top 
and  the  average  value  of  pressure  drop  during  the  purge  pro¬ 
cess  increases  with  the  operating  pressure.  Since  there  will  be 
the  second-falling  decreasing  stage  for  relatively  high  operating 
pressure,  like  at  2.00  bar,  the  purge  time  can  be  chosen  appro¬ 
priately  based  on  the  rules,  which  is  probably  deducted  from 
Fig.  1 1 ,  to  improve  the  purging  efficiency  by  prolonging  the  peri¬ 
ods  with  high  pressure  drop  before  the  second-falling  decreasing 
stage  comes. 

Fig.  12  shows  the  pressure  swing  distribution  along  the  chan¬ 
nel  during  the  process.  Similarly  with  the  pressure  drop  varia¬ 
tion,  the  pressure  swing  value  at  every  calculated  point  increases 
with  the  operating  pressure.  The  distribution  of  pressure  swing  at 
a  lower  operating  pressure  seems  a  little  homogeneous  than  the 
one  at  a  higher  operating  pressure.  However,  the  uniform  trend 
of  pressure  swing  variation  for  all  the  cases  is  that  the  value 


Fig.  12.  Pressure  swing  vs.  fractional  location  along  the  anode  channel  during 
the  purge  process  for  three  different  operating  pressures:  1 .50, 1 .75  and  2.00  bar. 
Aa  =  1  at  1.0  A  cm-2;  purge  time,  0.7  s. 
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increases  along  the  channel,  which  is  even  more  pronounced  at 
higher  operating  pressures.  Pressure  swing,  as  explained  above, 
reflects  the  highest  impulse  of  pressure  at  certain  point  and 
mechanical  effects  on  the  membrane  electrode  assembly.  There¬ 
fore,  the  area  of  the  MEA  close  to  the  exit  of  flow  field  always 
suffers  the  top  mechanical  attacks  of  gas  flow  in  this  transient 
process,  although  this  effect  seems  uniform  in  some  degree  in 
the  channel  at  low  operating  pressure. 

3.4.  Purge  time  effects 

The  length  of  purge  time  is  a  key  factor  to  be  considered  when 
the  purging  strategy  is  made.  Experiments  and  the  correspond¬ 
ing  simulations  with  the  same  operating  pressure  2.00  bar,  and 
current  density  1  A  cm-2,  but  with  different  lengths  of  purge 
time:  0.36,  0.68  and  1.0  s,  were  conducted  to  test  the  effect  of 
the  length  of  purge  time  on  pressure  drop  and  pressure  swing.  In 
Fig.  13,  it  can  be  seen  that  the  profile  of  pressure  drop  curve  is 
completely  shaped  by  the  length  of  purge  time.  The  whole  profile 
of  the  pressure  drop,  including  the  stable  and  the  second-falling 
periods,  can  be  observed  if  a  relatively  long  purge  time  (1.0  s) 
has  been  chosen,  while  if  a  short  purge  time  (0.36  s)  is  adopted, 
only  a  part  of  the  whole  profile  without  the  second-falling  period 
can  be  available.  Fig.  14  displays  the  comparison  of  the  pres¬ 
sure  swings  for  these  three  cases.  As  can  be  seen  from  the  figure, 
for  the  short  period  there  is  an  obvious  monotonic  increment  of 
pressure  swing  along  the  channel.  Although  the  distribution  of 
pressure  swing  in  the  whole  flow  field  is  non-uniform,  the  area 
upstream  of  MEA  receives  a  low  pressure  swing,  or  mechanical 
effect  of  the  gas  flow.  Oppositely,  for  the  long  purge  time  the 
distribution  of  pressure  swing  seems  more  uniform,  compared 
with  the  one  for  a  short  purge  time,  however  the  area  upstream 
of  MEA  may  also  suffer  a  higher  pressure  swing  which  means 
more  severe  attacks  of  gas  flow  to  the  MEA. 

As  indicated  above,  a  too  short  purge  period  will  lead  to  the 
non-uniformity  of  distribution  of  pressure  swing  or  mechanical 


Fig.  13.  Pressure  drop  curves  during  the  purge  process  for  three  different  lengths 
of  purge  time:  0.36,  0.68  and  1.0  s.  A.a  =  l  at  1.0  Acm-2;  operating  pressure, 
2.00  bar. 


Fig.  14.  Pressure  swing  vs.  fractional  location  along  the  anode  channel  during 
the  purge  process  for  three  different  lengths  of  purge  time:  0.36,  0.68  and  1.0  s. 
A.a  =  1  at  1.0  Acm-2;  operating  pressure,  2.00 bar. 


impulse  of  gas  flow  on  MEA,  and  the  high  pressure  drop  period 
cannot  be  made  full  use  of  to  gain  an  efficient  purge  process, 
while  a  too  long  purge  time  may  include  some  periods  with  low 
pressure  drop  that  would  be  useless  for  a  better  purge  process, 
and  the  area  upstream  of  MEA  has  to  suffer  high  pressure  swings 
and  severe  attacks  of  gas  flow,  which  may  damage  the  MEA  if 
they  take  place  frequently.  Therefore,  the  length  of  purge  time 
should  be  chosen  appropriately  for  an  efficient  purge  process 
with  low  potential  damaging  to  the  components. 

3.5.  Hydrogen  stoichiometric  ratio  effects 

All  the  simulations  and  experiments  discussed  above  were 
conducted  in  condition  that  the  hydrogen  outlet  was  closed 
down  and  the  purge  processes  were  accomplished  by  the  elec¬ 
tromagnetic  valve,  thus  the  hydrogen  stoichiometric  ratio  Aa  =  1 . 
Actually  in  the  application  of  fuel  cell,  the  hydrogen  stoichio¬ 
metric  ratio  is  not  always  1  when  the  hydrogen  gas  circulation 
is  adopted.  The  following  discussion  displays  the  comparison 
of  the  purge  processes  with  different  hydrogen  stoichiometric 
ratios. 

Experiments  with  the  same  operating  pressure  2.00  bar  and 
current  density  0.6  Acm-2,  but  with  different  outlet  boundary 
conditions  and  hydrogen  stoichiometric  ratios  (1,  1.1  and  1.2), 
were  conducted  to  test  the  effect  of  stoichiometric  ratio  on  the 
pressure  drop  and  the  pressure  swing  in  the  purge  processes.  It 
can  be  expected  that  when  the  stoichiometric  ratio  is  reset,  the 
whole  profile  of  the  pressure  variation  during  the  purge  process 
may  be  modified  because  the  hydrogen  gas  circulation  is  added 
and  the  outlet  boundary  condition  has  changed.  Therefore,  the 
validations  should  be  conducted  specially,  which  are  shown  in 
Fig.  15.  There  can  also  be  good  agreements  between  experiments 
and  simulations  after  some  model  modifications. 

In  Fig.  15,  it  can  be  seen  that  when  the  hydrogen  circulation 
is  added  and  hydrogen  stoichiometric  ratio  is  enhanced  from  1 
to  1.1  or  1.2,  the  second-falling  periods  have  disappeared  and 
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Fig.  15.  Comparison  of  experimental  and  model  pressure  curves  during  the 
purge  process  for  three  different  stoichiometric  ratios:  Aa  =  1,  1.1  and  1.2  (from 
left  to  right).  Current  density,  0.6  A  cm-2;  purge  time,  0.63  s;  operating  pressure, 
2.00  bar. 


Fig.  17.  Pressure  swing  vs.  fractional  location  along  the  anode  channel  during 
the  purge  process  for  three  different  stoichiometric  ratios:  Xa  =  1,  1.1  and  1.2. 
Current  density,  0.6  A  cm-2;  purge  time,  0.63  s;  operating  pressure,  2.00 bar. 


the  pressures  decrease  more  dramatically  during  the  purge  pro¬ 
cess.  Fig.  16  shows  the  comparison  of  pressure  drop  for  the  three 
cases.  From  this  figure,  it  can  be  seen  that  when  the  stoichio¬ 
metric  ratio  is  set  beyond  1,  the  pressure  drop  curves  display 
the  sharp  decreases  after  the  top  pressure  drop  value  other  than 
the  period  with  a  high  value  of  pressure  drop  before  the  second- 
falling  stage  as  the  situation  of  the  ratio  is  set  1.  In  Fig.  17, 
the  pressure  swing  curves  for  these  cases  are  given.  When  the 
stoichiometric  ratio  is  1 ,  there  is  a  monotonic  increment  of  the 
pressure  swing  along  the  flow  channel,  thus  the  area  upstream  of 
MEA  may  receive  a  relatively  low  mechanical  effect  of  gas  flow. 
Oppositely,  the  distribution  of  pressure  swing  in  the  flow  field 
become  more  uniform  when  the  stoichiometric  ratio  is  enhanced 
to  1 . 1  or  1 .2,  but  in  these  situations  the  area  upstream  of  the  MEA 


Fig.  16.  Pressure  drop  curves  during  the  purge  process  for  three  different  stoi¬ 
chiometric  ratios:  A.a=  1,  1.1  and  1.2.  Current  density,  0.6  A  cm-2;  purge  time, 
0.63  s;  operating  pressure,  2.00  bar. 


has  to  suffer  relatively  high  pressure  swing  and  severe  mechan¬ 
ical  effect  of  gas  flow. 

Therefore,  for  the  purge  process,  enhancement  of  the  stoi¬ 
chiometric  ratio  inevitably  removes  the  period  with  high  pres¬ 
sure  drop  before  the  second-falling  stage  that  can  be  used  to 
improve  the  purging  efficiency.  Although  higher  stoichiomet¬ 
ric  ratios  can  lead  to  more  uniform  distributions  of  the  pressure 
swing  and  the  mechanical  effects  of  gas  flow  on  the  whole  MEA, 
it  does  improve  the  potential  mechanical  attacks  of  gas  flow  to 
the  area  upstream  of  MEA,  compared  with  the  distribution  with 
stoichiometric  ratio  1 . 

3.6.  Current  density  effects 

For  the  fuel  cells  on  vehicles,  different  current  densities  are 
always  required  due  to  various  loading  requirements,  thus  the 
purge  process  may  be  conducted  under  varied  power  outputs. 
Current  density  reflects  the  consumption  of  hydrogen  gas,  which 
can  affect  the  pressure  in  the  flow  field.  The  following  simula¬ 
tions,  with  the  same  operating  pressure  2.00  bar,  stoichiometric 
ratio  1  and  purge  time  0.64  s,  but  with  different  current  densities: 
J- 0.6,  0.8  and  1.0  A  cm-2  were  conducted  to  test  the  effect  of 
current  density  on  pressure  drop  and  pressure  swing  during  the 
purge  process. 

Fig.  18  displays  the  pressure  drop  curves  under  different  cur¬ 
rent  densities.  In  the  comparison,  it  can  be  indicated  that  the 
operating  parameter — current  density,  has  no  apparent  effect 
on  the  pressure  variation  during  this  transient  process,  espe¬ 
cially  after  the  top  pressure  drop  value  and  the  curves  almost 
match  together.  This  is  probably  because  the  purge  process  is 
accomplished  in  a  very  short  period,  during  which  the  gas  con¬ 
sumption  is  too  low  to  show  an  obvious  effect  on  the  pressure 
in  the  channel,  compared  with  the  compensation  of  high  pres¬ 
sure  reactant  from  the  gas  source.  In  Fig.  19,  the  comparison  of 
pressure  swing  for  the  three  cases  is  given.  It  can  be  seen  that 
three  curves  have  similar  trends — the  pressure  swing  increases 
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Fig.  18.  Pressure  drop  curves  during  the  purge  process  for  three  different  cur¬ 
rent  densities:  J= 0.6,  0.8  and  1.0  A  cm-2.  ka  =  1;  purge  time,  0.64 s;  operating 
pressure,  2.00  bar. 


Fig.  19.  Pressure  swing  vs.  fractional  location  along  the  anode  channel  dur¬ 
ing  the  purge  process  for  three  different  current  densities:  /=  0.6,  0.8  and 
1.0  Acm-2.  ka  =  1;  purge  time,  0.64s;  operating  pressure,  2.00 bar. 

monotonically  along  the  channel.  The  curve  under  the  current 
density  of  J-  0.8  Acm-2  is  slightly  lower  than  others  just  due 
to  the  slightly  lower  initial  operating  pressure,  which  lowers  the 
whole  pressure  swing  curve  in  some  degree. 

As  indicted  above,  the  current  density  does  not  affect  the  pres¬ 
sure  variation  obviously  during  purge  processes.  Therefore,  for 
vehicle  use  purging  under  different  loading  conditions  makes  no 
difference  in  terms  of  improving  purging  efficiency  and  allevi¬ 
ating  the  mechanical  effects  of  gas  flow  or  potential  damages  to 
the  MEA. 


4.  Conclusions 

A  one-dimensional  computational  fluid  dynamics  model  of 
a  PEM  fuel  cell  was  developed.  The  transient  behavior  of  pres¬ 
sure  in  the  channels  during  purge  processes  has  been  studied. 
The  results  indicated  that  the  distributed  pressure,  pressure  drop 
and  pressure  swing  all  increased  with  the  increment  of  operating 
pressure.  With  a  high  operating  pressure,  the  pressure  drop  pro¬ 
file  showed  a  second-falling  stage  that  can  be  made  full  use  of  to 
improve  purging  efficiency,  but  the  pressure  swing  displayed  a 
homogeneous  distribution  with  a  relatively  low  operating  pres¬ 
sure.  A  long  purge  time  allowed  for  the  whole  part  of  the  pressure 
drop  profile,  while  a  short  purge  time  only  gave  a  part  of  the  pro¬ 
file  and  a  relatively  non-uniform  distribution  of  pressure  swing. 
When  the  stoichiometric  ratio  exceeded  1,  the  pressure  drop 
curve  changed  and  decreased  more  sharply  after  the  top  value, 
but  the  pressure  swing  displayed  a  more  uniform  distribution. 
Current  density  did  not  affect  the  pressure  drop  and  the  pressure 
swing  obviously  during  this  transient  process.  All  the  distribu¬ 
tion  rules  of  related  parameters  deducted  from  this  study  will  be 
helpful  for  optimizing  the  purging  strategies  on  vehicles. 
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